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vN Abstract. In this paper, we address the problem of orientation that natu- 

^ rally arises when representing shapes like curves or surfaces as currents. In the 

^1^ field of computational anatomy, the framework of currents has indeed proved 

-^^ very efficient to model a wide variety of shapes. However, in such approaches, 

orientation of shapes is a fundamental issue that can lead to several draw- 

Cn backs in treating certain kind of datasets. More specifically, problems occur 

Cn with structures like acute pikes because of canceling effects of currents or with 

data that consists in many disconnected pieces like fiber bundles for which 

' . ' currents require a consistent orientation of all pieces. As a promising alter- 

\^ native to currents, varifolds, introduced in the context of geometric measure 

r ^ theory by F. Almgren, allow the representation of any non-oriented manifold 

• (more generally any non-oriented rectifiable set). In particular, we explain how 

^ varifolds can encode numerically non-oriented objects both from the discrete 

I ^ I and continuous point of view. We show various ways to build a Hilbert space 

structure on the set of varifolds based on the theory of reproducing kernels. 

^_^ We show that, unlike the currents' setting, these metrics are consistent with 

^ shape volume (theorem |4.1| and we derive a formula for the variation of metric 

QQ with respect to the shape (theorem |4.2| ). Finally, we propose a generalization 

^^ to non-oriented shapes of registration algorithms in the context of Large De- 

^^ formations Metric Mapping (LDDMM), which we detail with a few examples 

\ r^ in the last part of the paper. 
• 
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1. Introduction 



Statistical shape analysis, the study of shape variability among a group of sub- 
jects, has been an active research field since the early works of [9], [17] and [19]. 
Important applications have been found notably in medical imaging leading to the 
^^ growing community of computational anatomy. The general approach of computa- 

Cd tional anatomy is to consider shape variability between subjects as resulting from 

deformations of the ambient space that one tries to estimate. Thus, a fundamental 
requirement for most algorithms is to have a relevant measure of residual shape 
difference. One of the main feature in this domain, though, is the wide variety of 
different shape structures coming from data measurements : images, landmarks, 
unlabeled distributions of points, curves, curve bundles, surfaces... These last ex- 
amples have been gathering many efforts, mostly due to the difficulty to design 
satisfying data attachment distances. In that regard, a particularly elegant frame- 
work has been proposed in ^15^, introducing in computational anatomy the notion 
of mathematical currents, which we briefly present in section [2J The construction 
of Hilbert metrics on currents' spaces has lead to robust distances between shapes 
that are also very convenient from a numerical point of view, which has proved 
of high interest in terms of application to shape registration between curves and 
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surfaces ([16],[IT]). 

Despite such important successes in computational anatomy, several critics have 
recently emerged related the use of currents for shape analysis. In our point of 
view, most of them are the consequence of the orientation dilemma that is totally 
inherent to the modelling of shapes as currents. We have gathered a synthesis of 
those problems in section [2j The central question would be to provide a repre- 
sentation that is completely free of shape orientation. This issue, however, is not 
new in the field of geometric measure theory from which the notion of currents 
originally appeared. Alternative concepts have been introduced and studied from 
a theoretical angle, notably the idea of varifolds suggested by F.Almgren and later 
developed by W.Allard. The purpose of this paper is precisely to adapt varifolds 
to the problematics of computational anatomy, which are fairly different from geo- 
metric measure theory. This is not the only possible approach to mention though, 
since the concept of normal cycles has also drawn recent attention in shape analysis 
[14,[8.. 

In the third section of the paper, we first present the mathematical object of var- 
ifolds as measures on the Grassmannian bundle, in the spirit of [2 . The real new 
challenge compared to geometric measure theory is to define metrics between vari- 
folds that are simultaneously 'fitted' to the comparison of non-oriented shapes and 
easily computable from a practical point of view. As for currents, a convenient way 
is given by the reproducing kernels' setting. Section [4] is therefore focused on the 
definition of kernels for the space of varifolds and shows the fundamental result of 



theorem 4.1 We also give a theoretical formula of variation of such metrics with 
respect to the geometrical support of shapes in theorem |4.2[ Finally, we present 
a first set of applications of our framework with an extension of classical large 
deformation matching algorithm to unoriented curves and surfaces. 

2. Representation of shapes with currents : strengths and 

limitations 

2.1. Currents in computational anatomy. In the following, we shall briefiy 
present the approach of currents in computational anatomy as first introduced in 
this field in [15] and developed later on in [10 . The essential features of currents, 
as we shall see, is that it provides an embedding of the set of all shapes of given 
dimension into a common Banach space and gives an intrinsic representation, in 
the sense that it is independent of shape parametrization. 
Throughout the text, we will adopt the following notations and definitions : 

• E : the embedding n-dimensional euclidean space of shapes. 

• A^E {0 < p < n) : p-times exterior product of E^ which is a vector space 
of dimension (^) spanned by the set of simple p- vectors ^i A .. A ^p. 

• A^£^ is equipped with the usual euclidean metric given for two simple p- 
vectors ^ = ^i A .. A ^p and r] = r]i A .. A rjp hy the determinant of the Gram 
matrix (^,77) = det{{£^i^r]j))ij. In particular, |^| gives the volume of the 
corresponding parallelotope. 

• ^o(^) •— Cq^E^K^E*) : the set of continuous p-dimensional differential 
forms on E vanishing at infinity, ^o(^) equipped with the infinite norm is 
thus a Banach space. 

This leads to the following definition of p-currents on E : 
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Definition 2.1. ^^{Ey , the space of all continuous linear forms on the space of 
differential forms, is by definition the space of p- currents on E. 

Note that in the particular case p = 0, the previous definition is exactly the one 
of usual distributions on E (dual of the space of real- valued functions). Just as 
for distributions, simplest examples of currents are given by Dirac currents 5^ 
with X G £^, ^ G hPE such that for any differential for uo G Vt^{E)^ we have 

Now, the relationship between shapes and currents lies fundamentally in the fact 
that every d-dimensional oriented sub-manifold X of E of finite volume can be 
represented by an element ofQQ{Ey. It is indeed a classical result from integration 
theory (|l3],[20]) that any (i- dimensional differential form can be integrated along 
X and thus : 



(1) Cx: UJ^ [ 



X 



defines an element of QJ^i^Ey . The space of d-currents is much larger though and 
contains many other interesting objects, among them rectifiable subsets of E^ which 
are the generalized submanifolds in the point of view of measure theory, can be 
treated as currents in the exact same way as submanifolds, by integrating differen- 
tial forms : such currents are called rectifiable currents in the literature (cf [13] and 
[22). 

Another major interest of the previous is the consistency of the currents represen- 
tation between the continuous and the discrete setting. Any discrete shape given as 
a set of points with a mesh can be transcribed into a current by associating to each 
simplex of the mesh a Dirac current located at the center with a simple d- vector en- 
coding the local volume element. It can be shown that this discrete representation 
converges to the continuous one given by equation [l] as the mesh becomes finer : 
we refer to [10] for more details. In numerical terms, this means that a shape can 
be encoded as a finite collection of point positions and (i- vectors. 
The last important element to address is the construction of a metric on the space 
of currents that can be both computed explicitly (at least for discrete surfaces) 
and that transcribes as well as possible the intuitive idea of closeness between two 
shapes. A particularly nice framework was proposed in [15 by joining currents 
with the theory of reproducing kernel Hilbert space (RKHS). This approach 
consists in defining a vector kernel on E {K : E x E ^ C{ApE)) and its associated 
RKHS W. It is then shown that, under some regularity assumptions on the kernel, 
the space of p-currents is continuously embedded in the dual W^ which is also a 
Hilbert space. As a result, currents and therefore shapes can be compared using 
a Hilbert norm which has the additional property that the inner product between 
two Dirac currents writes : 

allowing explicit computations when considering discrete shapes represented as fi- 
nite sums of diracs. Such norms can be also proved to have several interesting 
theoretical properties. For instance, kernel norms of a rectifiable current is domi- 
nated by its Hausdorff measure (or its d-dimensional volume) : 

Proposition 2.1. // Cx is a rectifiable current associated to a d— dimensional 
rectifiable set X, then : 

WCxWw < Cte-H^X) 
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where Cte is a constant depending only the choice of the kernel. 
Proof. For all uj eW, we have : 

ICxHI < / llculloo 

Jx 

< \cj\oon^{x) 

< Cte\\uj\\wn'^{X) 

where the last inequahty results from the continuous embedding of ^o(^) ^^^^ W. 
Thus, \\Cx\\w' <Cten^{X). D 

However, due to orientation, the converse domination does not hold, as we shall see 
in the following. 

Such kernel metrics on currents have been successfully used to derive data attach- 
ment distances, allowing an extension of Large Deformation Diffeomorphic Metric 
Mapping (LDDMM) algorithms for registration of curves or surfaces ( [H] , [16] , [TT] ) 
as well as statistical estimation of templates. Nevertheless, currents are intrinsi- 
cally modelling oriented objects : we shall show in the following subsection the 
limitations that this induces in terms of classical applications in shape analysis. 

2.2. Currents and the issue of shape orientation. As mentioned previously, 
the current associated to a given d— dimensional rectifiable subset of E depends on 
the orientation of X. If for instance X is connected and X denotes the same set but 
with opposite orientation then, in the space of currents, Cx = —Cx- It results that, 
in order to compare two shapes in the currents' setting, it is absolutely compulsory 
to have orientations of both shapes which must be in addition consistent with each 
other. A second point is that currents with opposite orientation within a small 
space domain may cancel each other with respect to the kernel metric. Indeed, for 
the case of a Gaussian kernel with scale a : 

(2) ||5«_5-«||^,=2|ep(l-e-l--^lV-^) 

which vanishes whenever \x — y\ is small compared to a. This trivial fact has im- 
portant consequences in practice. 

Difficulties can ffist occur with non-orient able shapes on which integration of differ- 
ential forms like in equation (fTj) is not well-defined and thus currents cannot handle 
unorientable shapes in a satisfying way. This ffist drawback remains anecdotal 
though, since most shapes in computational anatomy are orient able. However, 
even if orient able, orientating shapes consistently can be either a difficult or even 
ill-posed problem in certain datasets. This is notably the case for fiber bundles in 
the 2D or 3D space consisting of many different and possibly disconnected pieces 
of curves. Indeed, if a given shape has N connected components, there are 2^ dif- 
ferent possible orientations. In the example of figure [l] we have shown what could 
occur in terms of matching if bundles of small curves are oriented randomly. Using 
currents with such objects requires a way to propagate orientation from one part of 
the shape to all the rest. Technically, this leads to additional pre-processing which 
can be particularly tricky for highly disconnected bundle of curves with many dif- 
ferent directions like the ones corresponding to white matter fibers estimated from 
DTI in brain imaging, we show the example of one subject in figure [2] This example 
also emphasizes another difficulty appearing with sets of curves crossing each other 
for which the very notion of 'consistent' orientation may become meaningless. In 
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Figure 1 . Example of a matching between two fiber bundles con- 
sisting of many different pieces of curves, each of them being ori- 
ented randomly (left figure) versus a consistent orientation of all 
pieces (right figure). The source shape is in dark blue, the target 
in red and the deformed source is in light blue. We see that a 
random orientation does not provide a satisfying matching result 




Figure 2. An example of white matter fiber bundle estimated 
from Diffusion Tensor Imaging (DTI) illustrating the potential dif- 
ficulty of consistent orientation of all different fibers. 



such cases, the orientation is clearly irrelevant and one would like to treat objects 
as sets of unoriented shapes, which is impossible with currents. 

Finally, orientation may represent an obstacle even in the simplest case of usual 
oriented and connected submanifolds because structures like sharp spines or tails 
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naturally lead, when represented in a kernel Hilbert space of currents, to annihila- 
tion of some parts of the shapes. This comes again from equation ([2|, or, to state 
it in a more general way, from the fact that, while the I^-norm of a set is always 
dominated by its Hausdorff measure as we saw in proposition |2.1[ there exists sets 
of given Hausdorff measure but with arbitrarily small W^-noim.. The plane curve 
of figure [3] is an example. Indeed, the 1-Hausdorff measure being the length of the 
curve : ^^(7) = 2(1 + e) > 2. On the other hand, if 7 is parametrized on an interval 
/, then we have : 

l|C^II^'= / 7'{sfK{j{s),j{t)Mt)dsdt 

Jlxl 
Let's assume, to simplify computations, that kernel i^ is a translation invariant 
kernel of the form K{x^y) = k{x — y) Id^2 with Lipschitz regularity. Then, by 
expanding the integrals as the sum of the four pieces of curves 7^^ , 7^2 5 7^i 5 7^2 5 
we can eventually simplify the result thanks to the translation-invariance, which 
gives : 

Since \\C^^^ ||^, = \\C^^^ ||^, = //[o,i]x[o,i] Kl{^)-lif))^^dsdt, it is straightforward 
that ||C^,J|^, - {C^^^.C^Jw =0{e^). Moreover, 

\\C^., \\w'-{C^.,.C^.M' = II [^(7W - 7W) - ^(7(^) - lit) - (6,0))] dsdt 

^-'[0,l]x[0,l] 

Since the kernel k is assumed to be Lipschitz, 

Hlis) - 7(t)) - fc(7(s) - l{t) - (e,0))] < Ctee 

Therefore, we have eventually proved that ||C^||h^/ = 0{^/e) whereas for all e, 
^^(n) ^ 2. In practical terms, this means that some meaningful structures of 
shapes, modeled as elements of a RKHS of currents, can vanish in this representa- 
tion. We will detail, in a coming section, what important limitations this can induce 
in matching certain examples of curves or surfaces having such kind of structures. 

3. Varifolds 

In this section, we present varifolds first as a theoretical object, following the 
original work of Almgren in ^ and later of Allard in ^. We then explain in detail, 
in a spirit very similar to currents, the fundamental relationship between varifolds 
and unoriented rectifiable sets and show how this representation can be efficiently 
transcribed in a practical way. The construction of kernel metrics on varifolds is 
postponed to the next section. 

3.1. Grassmannian, varifolds : definitions and basic properties. Varifolds 
have been first introduced in the context of geometric measure theory as a con- 
venient way to address Plateau's problem of finding least area surfaces with a 
prescribed boundary. These developments clearly do not enter in the scope of this 
paper. Our purpose here is rather to connect conveniently varifolds to the frame- 
work and language of computational anatomy. Therefore, in the following, we will 
focus mainly on definitions of such objects and explain in what respect these defi- 
nitions are computationally relevant. The basic idea behind varifolds is to represent 
any rectifiable set (orientable or not) as a distribution of unoriented tangent spaces 
spread in the embedding space E. As we shall see, varifolds encompass not only 
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Figu re 3 . Counter-example to the reciprocal inequality of propo- 

For all e, ^^(7) = 2(1 + e) > 2 whereas \\C^\\w' — > 



sition 



2.1 



e^O 



manifolds and rectifiable sets but more generally sets of directions in the space. 
Let's take again the conventions of section [2J E being a vector space of dimension 
n and d an integer with < d < n. What we first need is a way to represent 
tangent spaces of dimension d in the space E. For this, we introduce briefly, in the 
following, Grassmann manifolds as it will prove useful later on. 

Definition 3.1. The Grassmann manifold of dimension d in E, denoted Gd{E), is 
the set of all d- dimensional subspaces of E. It can be identified to the quotient space 
of all families of d independent vectors of E by the equivalence relation obtained by 
identifying families that generate the same subspace. 

As a quotient space, one can show that Gd{E) inherits a structure of compact 
Riemannian manifold of dimension d{n—d) (cf [5, and [27^). Another way to think of 
an element in the Grassmann manifold is to identify a subspace V to the orthogonal 
projection on V, which embeds Gd{E) into the space of linear homomorphisms 
jC{E). For the following, we shall also need a practical way to think of the tangent 
spaces of Gd{E) and compute variations within this manifold. If V is an element 
of Gd{E), we can consider Uy = {W G Gd{E) | W n V^ = {0}}, which is an open 
subset of the Grassmann manifold. Every element in W G Uy can be written as 
W = {v ^l{v)^ '^ G V} for a certain linear function / G C{V^ V^) and / is uniquely 
determined by W^ as proved in [23|. Thus, we get a bijective map ijjv • l^v -^ 
£(V, V^). Referring again to [23], (^/y, ipy) defines an atlas on the manifold Gd{E) 
and, as a consequence, there is a natural isomorphism between the tangent space 
TyGd{E) and C{V^V^). Now, if for t g] — e,e[, V{t) is a differentiable curve on 
Gd{E) with V{0) = Vq^ how do we represent the derivative V'{0) ? Considering any 
curve v{t) with v{t) G V{t) for ah t, it is easy to see that V'{0).v{0) = py±{v'{0)) 
and this gives the definition of V^{0) as an element of jC{Vo^ ^o~)' 

Finally, let's mention a last way to embed the Grassmannian which is closer to 
the approach of section [2J This is formulated by the Pliicker embedding property 
below : 
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Proposition 3.1. The following application 

ip: Gd{E) -^ P{A^E) 

vect{vi,...,Vd) I — > [vl^...^Vd] 

where P {^K^E^ is the real projective space of K^E, is an embedding. It is even a 
homeomorphism in the cases d = 1 or d = n — 1. 

As a result, we can think of an unoriented (i- dimensional tangent space as an el- 
ement of P (A^£^), which is a space of dimension (^) — 1, much bigger in general 
than Gd{E). 

Following closely [2 , varifolds are now defined precisely as distributions of unori- 
ented tangent spaces in E : 

Definition 3.2. A d- dimensional varifold on E is a Borel finite measure (or dis- 
tribution) on the product space E x Gd{E), i.e an element of Cq{E x Gd{E))' . 

Note that this differs from the original definition of varifolds given by Almgren 
but his definition is equivalent to the previous one as explained in the preface of [3] . 
Diracs in the space of varifolds are of the form Sf^^y) with x e E and V G Gd{E) 
and act on functions of Co{E x Gd{E)) by the relation : 

Vc^ e Co{E X Gd{E)), 5^,y^{cj)=cj{x,V) 

In this context, a dirac consists in the data of a position x in the space attached to 
a d-dimensional (non-oriented) subspace V that will play the role of tangent space 
in the case of rectifiable sets. 

3.2. Non-oriented shapes as varifolds. In section |2J we have seen that ori- 
ented rectifiable sets of dimension d are canonically represented as d-currents. In 
the case of non-oriented shapes, the right notion is precisely the one of varifolds. 
Indeed, let X be a non-oriented rectifiable set of E of dimension d. To X, one 
can associate a varifold /ix, which, in the measure point of view, is given by : 
l^x{A) = n'^iixeX \ {x.T^X) e A}) for ah Borel subset A C E x Gd{E). By 
Riesz representation theorem, jj^x has its unique equivalent in terms of continuous 
linear form on Co{E x Gd{E)). For ah ou e Co{E x Gd{E)) : 

(3) /ixM= / uj{x,V)d^x{x,V) 

JExGdiE) 

Such special varifolds are called rectifiable varifolds. For more explicit expression, 
assume that X is a smooth compact submanifold of E with a parametrization 
J :U ^ E with U an open subset of R"^, then for ah uj e Co{E x Gd{E)) 



(4) /ix(^) = / uj{j{u),r]{u))\j\u)\du 

Ju 

with the notations j\u) = Af^^d-f/dui G A'^(^), r]{u) = [-/'{u)] e Gd{E). Note 
that the integral of equation Q is, as expected, independent of any reparametriza- 
tion of X, positively or negatively oriented. Now, in discrete geometry, shapes are 
given as polygonal sets of points which are also encompassed in the category of 
rectifiable subsets and therefore represent able as varifolds. In the same spirit as 
section |2J any mesh set can be coded as a finite sum of dirac varifolds of the form 
^7=1 ^i-^iPi^Vi)- Again, n is the number of cells of the mesh, pi G E are the centers 
of each cell, Vi G Gd{E) the unoriented tangent space to the shape at point Pi 
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and 



Ci e 



the volume of the i-th ceh. More specificahy, let's examine the two 



most usual examples in practice. First, consider a curve X given as a set of points 
{xk}k=i..N and a set of edges {fl^ff) G {1, .., A/"}^ for i = 1, ..,n. Then, with the 
previous conventions : 



(5) 



Pi 



Xfl + Xf2 



Vi = {xfiXfi) 



XflXf2 



frf? 



Note that all the previous equations remain unchanged if for any i we take (/^^, //) 
instead of (fl^ff) as a face, which is consistent with the idea of non-orientation. 
In a similar way, if X is now for instance a triangulated surface in £^ = R^ given 
by a set of points {xk}k=i..N with a set of n triangles (//, ff^ ff) G {1, .., N}^ for 
i = 1, .., n then we get : 



(6) 



Here 



Pi 



V,. 



Xfl + Xf2 + Xf3 
J i J i J i 



XflXf2^ XflXf3 
J i J % J i J i 



-.\\XfiXp A xpXf?>\ 



XfiXf2^XfiXf3 denotes the 2— dimensional space generated by Xfixp and 

XfiXfL Again, all these equations are not dependent on the orientation given to 
each triangle. 

A last important point for the following is to be able to express the transport 
of varifolds by diffeomorphism in a way that is compatible with the transport of 
shapes. Let's fix a varifold /i G Co{E x Gd^E))' and (j) G Diff(£^), we define the 
transport ^*/i of /i by ^ by pull-back and push- forward operations : 



(7) 



Vw G Co{E X Gd{E)), (</.,m) {<^) = M {<P*oj) 



where <p*uj is the pull-back of w by (p. For any x E E and V G Gd{E) a d-dimensional 
subspace with orthonormal basis {ui, ...,Ud), 4>*oj is defined by the relation : 



(8) 



*CJ)(X,V) = \d:,(l){ui) ^...^d:,(|){Ud)\ Uj{(j){x),d:,(j).V). 



In the last equation, dx4>y denotes the element of Gd{E) that is the image of V 
by dxCJ)^ the term \dx(l){ui) A ... A dx(l){ud)\ is the d-dimensional Jacobian of (j) on 
the subspace V and represents the local change of d-dimensional volume of the 
transformation. With equations ^ and ([8|, it is straightforward to compute the 
transport of a dirac St^^y^ : 



(9) 



(t>Jix,v) = \dx(l){ui) A ... A dx(l){ud)\ •5{(j,{x)4^(j).v) 



if (lii, ...^Ud) is an orthonormal basis of V. Now, with the previous definitions, we 
can show that the varifold representation of a rectifiable set commutes with the 
transport by diffeomorphism in the sense given by the following proposition : 
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Proposition 3.2. If X is a d- dimensional rectifiable subset of E and (j) G Diff{E), 
then : 

Proof. If X is a rectifiable subset of E^ there exists (cf [2^) at almost every point 
of X a tangent space TxX and we have for all uj : 

/ixM= / uj{x,TxX)dn\x) 
Jx 

and thus 



j {c^^uj){x,TxX)dn\x) 

JX 

uj{c^{x),dxC^(TxX))\dxC^.TxX\dn\x) 



ix 

where \dx(l).TxX\ is the (i- dimensional Jacobian of (j) along the tangent space T^X^ 
as previously. For almost all x G X, dx(l){TxX) is a tangent space of (/>(X) and by 
the generalization of the change of variables for rectifiable subsets (corollary 3.2.20 
in [l3]), we obtain : 






Thus, for any cj, (j)^iix{^) = M0(x)(^) and we have proved proposition 



3.2 



D 



In conclusion to this section, we have seen that varifolds offer a satisfying frame- 
work to represent non-oriented shapes both from the continuous setting and in the 
computational cases of meshed curves, surfaces... Up to this point, we have shown 
how to represent computationally a meshed shape as a finite set of diracs, each of 
them carrying a local information of position, unoriented tangent space and local 
volume and derive the equations of varifolds' transport by deformation. In the pur- 
pose of adapting large deformation matching to varifolds (section[5|, an appropriate 
metric still needs to be defined, a topic that is discussed in the next section. 

4. Kernel metrics on varifolds 

The use of reproducing kernels is a fundamental step when working with currents 
because it provides regularized metrics for the comparison of shapes which has the 
additional advantage of having an underlying Hilbert space structure. We have also 
mentioned in section |2] that kernels are particularly well-fitted to compute distances 
between discretized shapes because of the simple expression of dot product between 
two diracs. All these reasons motivate a similar approach in the treatment of 
unoriented shapes through the varifold setting that has been presented. It's worth 
mentioning that other approaches could be possible, by working with more general 
Riemannian metrics on distribution of tangent spaces, generalizing for instance 
what is done for the Grassmann manifold in pLj. We argue however that RKHS 
metrics are very convenient from a computational point of view, especially in our 
applications to shape matching and analysis. In this section, we propose a generic 
and simple way to build relevant reproducing kernels on the space of varifolds, by 
basically making tensor products of kernels on E and on Gd{E). 
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4.1. Kernels on the Grassmann manifold. We start by the construction of 
kernels on Gd{E). This problem is not totally new since examples of such kernels 
have been provided in several recent works related to machine learning, notably in 
[.18 and [26 . Such constructions are usually based on the notion of principle angles 
between two subspaces. An alternative way, which has the advantage of offering 
a very wide class of possible kernels on the Grassmannian, is to simply use the 
embedding of Gd{E) into C{E) mentioned in section [sj It identifies any subspace 
V G Gd{E) with the orthogonal projection on F, which we denote by pv- As a 
consequence, one can induce straightforwardly a kernel on Gd{E) by restriction of 
a positive kernel defined on C{E). Since jC{E) is a finite-dimensional vector space, 
there are no difficulties in defining such kernels. Typical choices are given by : 

(10) kcB{V,W) = {pv,Pw)'' 

with /c G N* and (., .) the usual Frobenius metric in C{E). More generally, any 
function P{{pv^Pw))^ P being a polynomial with positive coefficients, is a positive 
kernel on C{E). Other important possibilities are the kernels induced by Gaussians 
in jC{E), namely : 

_ \pv -Pw\'^ 

(11) kG{V,W) = e~ Cr2 

Such kernels allow the comparison of subspaces with respect to a certain scale given 
by parameter a. Many other could be defined with this method. It's also possible to 
express these kernels using the principle angles. If V and W are two d-dimensional 
subspaces of E^ the d principle angles ^i, .., ^^ between them are defined recursively 
by the relations : 

cosiOk) = max max (^ — ;,- r) wztri 

vkevwkew'\vk\ \wk\ 

Vz G {l,..,/c - 1}, Vk-Lvi, Wk-Lwi 

Thus, there exists orthonormal frames ('Ui, ..,Vd) and {wi^ --^Wd) of V and W such 
that cos{Oi) = {vi^Wi). As easily seen, see e.g. [2^, there are at most min((i, n — 
d) non-zero principle angles and they completely determine the relative position 
between the two subspaces. Now it's an easy verification that : 

d d 

(py,PH^) = 2^(^„^i)' = 2^cos2(^,) 

i=l i=l 



It results that, up to a scaling factor, the kernel of equation (10) has the following 
expression : 

(12) kcB{V,W)= (^cos2(^0) 

For /c = 1, we see that the kernel we obtain is precisely the classical Cauchy-Binet 
kernel on the Grassmannian, as in [^. The Gaussian kernel also has the alternative 
expression : 

(13) kG{V,W) = e-^-Sti(i-cos^(^.)) = J]e-^^"'(^*) 
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This formulation with principle angles offers, in our sense, a nicely interpretable 
way to understand the behavior of such kernels with respect to the subspaces' 
relative position. The general construction we have proposed thus provides a wide 
variety of induced kernels on Gd{E) that are effectively computable (either by the 
expression of the projection matrix or by computing the principle angles). 

Remark 4.1. We also want to mention, without entering into details, another 
possible way of building kernels on Gd{E), relying on the Pliicker embedding of 
section \3.1\ instead of the identification with projectors. Indeed, one can induce as 
previously a kernel on the Grassmann manifold from a kernel on P{A^E). Since 
this projective space can be also identified to the unit sphere of A^{E) quotiented 
by the reorientation group {±1}; it is not difficult to provide kernels on P{A^E) 
by defining kernels on the vector space A^E, invariant with respect to the action 
of {±1}. For curves or hypersurfaces, the Pliicker embedding being a homeomor- 
phism, this approach is relevant and provides kernels quite similar to our previous 
method. However, for high dimension and codimension, we see that the embedding 
space becomes very high- dimensional (^) — 1 compared to the actual dimension of 
Gd{E) and we believe that this can be a limiting factor both from a theoretical and 
computational point of view. 

4.2. Construction of kernels on varifolds. The issue of kernels on Gd{E) being 
addressed, we now move to the case of varifolds. Since we are considering functions 
defined on a product space, a natural way to build appropriate kernels is the tensor 
product trick. We first remind a well-known general property from reproducing 
kernels' theory: 

Lemma 4.1. Let A and B be two sets and kA, ks two positive real kernels respec- 
tively on A and B. Then, the tensor product kA ^ ks defined by : 

kA ^ ks ((a, 6), {a\ h')) = kA{a, a') ksib, b') 

is a positive real kernel on Ax B. 

Proof. For any A/" G N*, ai,..,aAr G A and 6i,...,6Ar G B then, since kA and ks 
are positive kernels, the matrices [kA{cii-,cij)]^ ■ and [^^(^i^^j)]^ • are positive by 
definition. The matrix [/c^ /^^((^^i, ^i), («j, ^j))]^ • being the Hadamard product 
of the two previous matrices, it follows from Schur product theorem that it is also 
positive. Thus, kA /c^ is a positive kernel. D 

We remind that to the kernel kA ^ ks is associated a unique Hilbert space of 
real- valued functions on Ax B (the RKHS of the kernel) on which all linear forms 
^(a,6) • / ^ /(<^7 b) are continuous. Now, going back to the case of varifolds, the 
following holds : 

Proposition 4.1. Assume that we are given a positive real kernel k^ on the space E 
such that ke is continuous, bounded and for all x G E, the function ke{x^ .) vanishes 
at infinity. Assume that a second kernel kt is defined on the manifold Gd{E) and is 
also continuous. Then the RKHS W associated to the kernel ke h is continuously 
embedded into the space Cq{E x Gd{E)). 

Proof. By definition : 

(14) k ((x, V), {y, V)) = Ux, y) hiV, V) 
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and so, thanks to the assumptions on the kernels, A:((x,V),.) is continuous on 
E X Gd{E) and belongs to Co{E x Gd{E)). W being the Hilbert space generated 
by the functions k{{x, V), .), we have W C Co{E x Gd{E)). Moreover, if uj e W : 

uj{x,V) =(5(^,y)(cj) = {k{{x,V),.),uj)w 

With Cauchy-Schwarz inequality : |a;(x,F)| < ||/c((x, V), .)||vF-||^||Ty- In addition, 
||/c((x, V), .)||h^ = ^/k{{x^V)^{x^V)) and both kernels ke and kt are bounded so 
that k is also bounded. We conclude that |a;|oo < vT^Uo-II^IIvf? which proves that 
the inclusion embedding i : W ^^ Co{E x Gd{E)) is indeed continuous. D 

Consequently, there exists a continuous mapping of the space of varifolds Co{E x 
Gd{E)y into the dual of W. Just as for currents, we can compare varifolds and 
unoriented subsets through the Hilbert norm of W . As a result of the kernel 
property, if xi,X2 G E and V\^V2 ^ Gd{E) : 

(15) (^(a^i,yi),^(x2,y2))H^' =^e(3:i,:r2)A:t(Vi,F2) 

We have therefore a generic way to build kernels for varifolds which are separable 
since such kernels are tensor products of a kernel on the ambient space E with a 
kernel on the set of all tangent spaces Gd{E). Building kernels on the euclidean 
space E raises no difficulties (examples were already given in section |2| and we 
have presented a way to build kernels on Gd{E). However, it still remains unclear 
whether the dual application i* : Co{E x Gd{E)y -^ W is always an embedding. 
In general, this is actually not the case : although W ^^ Co{E x Gd{E)) is an 
embedding, the dual application need not be injective. For instance, choosing the 



Cauchy-Binet kernel of equation (10) on the Grassmannian makes i* not injective. 
The fact that i* is an embedding is called the Co-universality property of the 
kernel and, as proved in [6 , it is equivalent to the property of W being dense in 
Co{E X Gd{E)). With respect to varifolds' kernels, we can still show a general but 
weaker result : under a few assumptions, z* is injective on the set of finite unions 
of submanifolds. The result is the following : 

Proposition 4.2. Let k = ke ^ kt be a kernel as in proposition \4-l\ Assume that 
kernel ke is Co-universal and that the kernel kt is such that kt{V^V) > for all 
V G Gd{E). Let X = Ui^i ^i ^^^ ^ — Uj^i ^j ^^ ^^^ finite union of compact 
d- dimensional submanifolds of E. If \\jiix — My||vF' = then X = Y. 

Proof We will denote by We and Wt the RKHS associated to kernels ke and kt. 
Let's start by the case where X is a single submanifold of £^. If \\/j^x—/J^y\\w' = and 
X 7^ y, one can find xq G X\Y and there exists r > such that 5(xo, r) fl F = 0. 
Let a G Co(£^,M) be any continuous function such that the support of a is included 

in B{xo^r). Let ujt = ktiT^^X^ .) G Wt. For all functions uje G We, cjg cjt G VF 
and thus (/ix — My)('^e ^ ^t) = 0, i.e 

/ Ue{x)ut{T,X)dn\x) - [ uje{y)ujt(TyY)dn\y) = 
Jx Jy 

Now, since We is dense in Co(£^,M), we can approximate uniformly a by functions 
in We. By uniform convergence in the previous integral and the fact that a = in 
y, we get : 

a{x)uJt{T^X)dn'^{x) = 



L 
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This holds for ah continuous functions a supported in 5(xo,r), which is clearly 
impossible since x \-^ uJt{TxX) is continuous and uotiTx^X) = ktiTx^X^Tx^X) ^ 0. 
If X is now a finite reunion of compact submanifolds and X ^Y then, as previously, 
we can find xq G Ui=i ^i with B{xo,ro) fl F = 0. B{xo^ro) fl X is itself a non- 
empty finite reunion of submanifolds. In order to get a similar proof as for the 
one submanifold case, we need to show the following : there exists a point x G 
B{xo^ro) n X, p > with B{x^p) C B{xo^ro) and j G {1,..,A^} such that for all 
function uj supported in B{x^p) x Gd{E)^ Px{^) = Mx (<^)- This can be proved 
recursively. 

Assuming the result for X — 1, we can find xi, ri,ji with B{xi^ ri) C B{xo^ tq) and 
/i. jAT-i^. = fiXj^ for functions supported in B{xi^ri). Now, there are two distinct 
cases : either B{xi,ri) fl Xj^ = B{xi^ri) fl Xn or we can assume for instance 
that there exists X2 G B{xi,ri) HX^^, r2 > such that B{x2^r2) C B{xi^ri) and 
B{x2^r2) n Xn = 0. In the first case, since Xn coincides with Xj^ on 5(xi,ri), 
we have still /iyAr j^, = juxj for functions supported in B{xi^ri) x Gd{E) and we 
take X = xi, p = ri, j = ji. In the second case, one can take ^ = ^2, P = ^2, j = ji 
and the result holds as well. The rest of the proof is then exactly similar to the 
case of a unique submanifold. D 



Even though there is no general injectivity for the dual application, proposition |4.2| 
ensures that the metrics we use are at least able to distinguish finite reunion of 
submanifolds. We believe that the result could be generalized to rectifiable subsets 
of E but the proof might be more involved and technical. 

In the specific case of the Gaussian kernel on Gd{E) given in section 4.1, we can 
actually recover the injectivity of z* on Co{E x Gd{E)y itself, which is stated in 
the next proposition : 

Proposition 4.3. If k = he ^ kc with ke a Co-universal kernel on E and kc 



the restriction of a Gaussian kernel on C{E) given in equation (11), then k is a 
Co-universal kernel on E x Gd{E). 

Proof The proof is mainly based on two results about RKHS. The first one is 
the very well-known property that a Gaussian kernel on a finite dimensional vector 
space is Co-universal. In our context, the Gaussian kernel K defined on C{E) x C{E) 

\ll-l2\'^ 

by K(li^l2) = e ^^2 is thus Co-universal. Now, the kernel kc that we have 



defined in equation (11) is the restriction of K to the subset of all orthogonal 
projections on c^-dimensional subspaces (which is identified to Gd{E)). This subset 
is closed in jC{E) and it is proved in [6 (corollary 3) that kc is then a Co-universal 
kernel. If We C Co(£^,M) and Wt C Co(Grf(£'),R) are the RKHS corresponding 
to kernels ke and /cg, it results that we have We dense in Co(^,M) (because ke is 
assumed to be Co-universal) and Wt dense in Co{Gd{E)^R). It is then clear that 
W, which is the completion of We ^ Wt, is dense in Co{E x Gd{E),R). D 

In summary, we have explained in this subsection how to define Hilbert metrics on 
varifolds that are computed from tensor products between kernels on E and kernels 
on the Grassmann manifold. This does not provide necessarily real distances on 
varifolds because of the possible non- injectivity of the application between the space 
of varifolds and the dual of the RKHS W. Yet, we have shown that for a very 
wide class of such kernels, the resulting distances are separating finite unions of 
submanifolds and we have argued in favor of Gaussian kernels on the Grassmann 



VARIFOLDS 15 

manifold for which we obtain the real injectivity and thus distances on the space 
of all varifolds. We will elaborate a little more on the interest of such Gaussian 
kernels in specific situations in the next subsection (cf figure pi. 

4.3. Properties of RKHS norms on varifolds. There are several interesting 
remarks and properties that we can state about the previously defined norms, no- 
tably if we compare it to the framework of currents that was exposed in the first 
section. First of all, let's consider the limit case of an infinite scale for the Gaussian 
kernel on the space E^ i.e ke{x^y) = Me- If X is any (i- dimensional submanifold 
of E^ the norm of the current Cx becomes : 

11^x11^. = ff {^ix),m) = { I iix). I i{x)) 

JJXxX Jx Jx 

Now, from Stokes' formula, J^ ^{x) is a term that only depends on the boundary 
of X. Therefore, at large scale, RKHS norms on currents are only sensitive to the 
boundaries of objects, and all shapes with no boundary vanishes in this represen- 
tation. 

In the varifold case, the behavior is fundamentally different. At large scale for 
the kernel /Ce, the RKHS distance becomes a distance between the distributions of 
tangent space directions depending on the kernel kt. More specifically : 

Wf^xfw = II kt{T,X,TyX)dn^{x)dn^{y) 

J JxxX 

= // kt{u,v)dux{u)dux{v) 

J JGd{E)xGd{E) 

where vx '-= H^ o Tan~^ is the image measure of 1-L^ by the application Tan : 
X -^ Gd{E), X \-^ TxX. In other words, at infinite scale for /cg, ||/ix||v^/ represents 
a metric on the distribution of the non-oriented tangent spaces to the shape X, 
which we can see, in a certain way, as a histogram metric on Gd{E). Therefore, 
even at infinite scale for ke , something of the shape still remains in the varifold rep- 
resentation contrarily to currents. In those situations, the choice of the kernel on 
the Grassmann part of varifolds is quite decisive. In particular, the Co-universality 
issue discussed previously becomes important. We show an illustration of such 
phenomenon in figure [4] for the two kernels explicitly given in section 4.1 : the 



Cauchy-Binet kernel of equation (10), which is not Co-universal and the Gaussian 
kernel of equation ( pT] ), which was proved to be. The conclusion that we can be 
drawn in general is that Gaussian kernels are obviously more effective in situations 
of multiple directions crossing at nearby points, as might appear when treating 
fibers or tree-Hke structures. 



We now come the problem raised by our discussion related to figure [3J We show in 
which precise sense RKHS metrics on varifolds successfully avoid the cancellation 
phenomenon of currents. This is formulated in the theorem below and its corollary : 

Theorem 4.1. Let ji he a positive d- dimensional varifold of E such that ja G W 
and Supp{ii) C 5(0, 1) x Gd{E). Let k be a reproducing kernel on E x Gd{E) and 



W its RKHS, defined like in section 4-2 by k{{x,u)^ iVi^)) = ke{x^y).kt{u^v) where 
ke is a reproducing kernel on E and kt a reproducing kernel on the Grassmannian. 
We make the following additional assumptions : 
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Figure 4. Computation of varifolds' distance between a 2D grid- 
like shape (in black) and its rotated version (in magenta) for vari- 
ous rotation angles and for both the Cauchy-Binet kernel and the 
Gaussian kernel on the Grassmann manifold, all in the case of a 
very large-scale kernel on the space E. The graph displays the ratio 
between the varifold distances and the norm of the original shape 
as a function of the rotation angle. With the Cauchy-Binet kernel, 
the values remain very small for all angles, and it is nearly unable 
to distinguish both shapes whereas the Gaussian kernel shows an 
expected behavior with a maximum distance for a 45° angle. 



1. ke is a radial scalar kernel so that we can write ke{x,y) = he{\x — y\). he is 
assumed to be a continuous positive function, with he{0) > 0. 

2. kt is a continuous and positive function on Gd{E) x Gd{E) such that for all 
ueGd{E), kt{u,u) >0. 

Then there exists a constant Cte independent of ji such that : 

M\w' > Cte^{ExGd{E)). 
Proof For any varifold /i G W\ we have : 

(16) iw;..^// ..(..ri.M«..wx,.)«,.(...)^ 

JJiExGdiE)^ 

For the proof of theorem |4.1[ we shall first examine the case of a constant kernel 
for ke. 

Step 1: We first assume that ke{x^y) = 1 for all x^y. Let's denote hy p : E x 
Gd{E) -^ Gd{E) the application (x^u) \-^ u. We introduce the image measure 
V '.= ji o p~^ defined on Gd{E). Then : 



Mfw = 11^ 



kt{u^v)du{u)du{v) . 

/GdiE)xGdiE) 

Note that v{Gd{E)) = fi{p~^{Gd{E))) = fi{E x GdE). Now, the compact group 
SO{E) of direct isometrics of E acts transitively on Gd{E) by the relation 
^.Span(ei, ..,6^^) := Span(^(ei), ..,^(6^^)), making Gd{E) a homogeneous space. 
Since SO{E) is a compact group, we can consider its unique bi-invariant Haar 
measure which we denote A with the convention X{SO{E)) = 1. If we fix a par- 
ticular element uq G Gd{E)^ A induces in turn a measure on Gd{E) defined for 
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all B by XcriB) = X{{g G SO{E) \ q.uq G B}). One can check that the right- 
invariance of A implies that Xqt does not depend on the choice of uq G Gd{E). In 
the same way, thanks to the left-invariance of A, Xcr is invariant by the action of 
G i.e Xcrig-B) = XGr{B) for all g. 

A second important element is that, as a homogeneous space under the action 
of a compact Lie group (on which there exists a bi-invariant metric), Gd{E) can be 
equipped by projection with a left- invariant distance with respect to the action of 
SO{E), which we will denote d. Let's now consider a continuous, positive L^ (for 
the measure Xcr) function cj) on Gd{E). Since XGr{Gd{E)) = X{SO{E)) = 1 < +00, 
(/) is also in L^ and we assume that \\(j)\\2 = 1- We can assume in addition that the 
support Supp(^) of (j) is included in a certain ball of radius S centered at ^^o- We 
also introduce the regularization function t/j on Gd{E) x Gd{E) defined by : 



tlj{u,v)= / (l){g.u)(t){g.v)dX{g). 

Jso(E) 



ISO{E 

It's obvious that '^ is also a positive continuous function and Cauchy-Schwarz in- 
equality gives : 

V(ii, v) e Gd{E)\ ^{u, v)< [ ^{g.ufdX{g) . 

JSO{E) 

The definition and invariance properties of A^^ show that for all u G Gd{E), if 
Tu '- g ^^ g-u then, A o r~^ = Xcr and it results that : 

/ ^{g.ufdX{g) = I ^{vfdXcM = UWl = 1 , 

JSO{E) JGd{E) 

therefore il){u^v) < 1. In addition, if d{u^v) > 2(5, we have d{g.u^g.v) > 2S for all 
g since d is left- invariant. Therefore d{g.u^uo) -\- d{g.v^uo) > d{g.u^g.v) > 26. It 
results that either d{g.u^uo) > 6 or d{g.v^uo) > S and anyhow (j){g.u)(j){g.v) = 0. 
Consequently, ip vanishes for d{u^ v) > 26. Now, since kt is positive 



M'w > 11^ 



kt{u^ v)iIj{u^ v)dh'{u)dh'{v) . 

Gd{E)xGd{E) 

Moreover, kt is also continuous and kt{u^u) > for all u so, by compactness of 
Gd{E), there exists a > and an open domain D C Gd{E) x Gd{E) of the form 
D = {{u^v) I d{u^v) < e} such that for all u^v G D^ ktiu^v) > a. We can 
also assume, by choosing appropriately the function ^, that S < ^ and therefore 
il){u^v) = outside D. It results that : 

(17) ||/i||^, ^ ^ // ^K v)dv{u)dv{v) . 

J JGd{E)xGd{E) 



18 NICOLAS CHARON, ALAIN TROUVE 

Using Fubini's theorem in the last integral, we can write 
// ilj{u^v)dh'{u)dh'{v) 



/Gd{E)xGd{E) 



[ i II ^{g.u)^{g.v)dv{u)dv{v) ) d\{g) 

JSO{E) \J JGd{E)xGd{E) J 

I (I ^{g.u)dv{u)\ dX{g) 

JSO{E) \JGd{E) J 

/ \ 2 



(18) > / / c^{g.u)dv{u)d\{g) 

\JSO{E) JGd{E) J 

the last estimate resulting from Cauchy-Schwartz inequality and the fact that 
X{SO{E)) = 1. Using again Fubini's theorem, 



/ / <l,{g.u)du{u)dX{g) = / ( / 

JsOiE) JGd{E) JGd{E) \JS 



(t){g.u)d\{g) dv(u) . 

SOiE) J 

Making the change of variable Tu - g ^-^ g.u in the inside integral and using the 
same argument as previously, we obtain that 

/ ^{g.u)dX{g) = I ^{w)dXGr{w) = Uh > 

JsOiE) JGd{E) 

Inserting the previous in equation (18), 



// 

J JG 



^{u,v)du{u)du{v) > \mlu{Gd{E))^ = \mlKE X GdEf . 

GdiE)xGdiE) 



So that we eventually obtain : 

Wtixfw, > a\\4>\\lfi{E X GdE)^ 

and since /3 := Q^||(/>||i is a constant that does not depend on /i, this concludes the 
proof in that case. 

Step 2: We now move to the proof for a general kernel /cg. From the hypotheses, 
the support of /i is included in 5(0, 1) x Gd{E) and ke is a continuous and positive 
radial scalar kernel on E. Thus, there exists a real number 5 < 1 such that for any 
x^y G E with \x — y\ < S^ we have k^ix^y) > he{0)/2 := k,. Let's now consider a 
covering of the unit ball with a set of cubes {Ci}i=i^,,^M of diameter smaller than 
S. Then the subsets Si := {Ci x Gd{E)) fi Supp(/i) form a partition of Supp(/i). 
Moreover, 



= // ke{x,y)kt{u,v)dii{x,u)dii{y,v) 

J J{ExGd{E)y 
M 

-yZ ^e(^' ^) ^^(^' ^)<^m(^' u)diii{y, v) 

i=l J JSiXS, 

M 

> n2_. II kt{u^v)djj.{x^u)djj.{y^v) . 
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We can apply the result of step 1 to each of these integrals, thus obtaining 

M 

the second inequality being a discrete Cauchy- Schwartz and the last one resulting 
from the fact that {Si} is a partition of Supp(/i). This ends the proof because the 
constants j3 and n are both independent of /i, and so is M which only depends on 
the kernel kp. D 



The hypotheses of theorem 4.1 on the kernels are mostly technical but are not 



particularly restrictive in practice since one can check easily that all the kernels of 



equations (10) and (11) comply to the requirements on kt. As a direct corollary, 



we also have the following result : 

Corollary 4.1. Let X he a rectifiable subset of E included in the unit ball. We 



make the same hypotheses on the kernel as in theorem \4.1\ Then there exists a 
constant Cte independent of X such that 

\\^x\\w^>Cten^{X). 



Proof This is essentially a special case of theorem |4.1| for rectifiable varifolds. In- 
deed, the support of fix is included in 5(0, 1) x Gd{E) and thanks to the theo- 
rem, ll/ixllvF' ^ Ctejiix{E X Gd{E)). In addition, from the very definition of jux^ 

fix{ExGd{E))=n^{X). D 

This result is theoretically essential because it shows that pathological cases as 
the one of figure [3] cannot happen when shapes are represented as varifolds. In 
practical applications, it ensures the consistency of the kernel norm we use with the 
actual volume of the shapes so that artificial elimination of mass during registration 
process or template estimation are not likely to occur in this setting. Note that 
when X is of codimension one and is defined as the boundary of a bounded domain, 
the notion of volume of a shape X is the area of the boundary of the domain and 
not the volume of the domain. In particular, the corollary says that when a plain 
shape is represented by its boundary, its norm for generic varifold kernel is larger 
than the area of its boundary. We will show more specific examples of this in the 
next section. 

A last fundamental issue to mention is the question of the variation of kernel 
norms with respect to the geometrical support of shapes, since registration and tem- 
plate estimation algorithms rely on the computations of gradients of such norms. 
We shall detail the technical computations of such gradients for discrete shapes in 
appendix B but it is also very valuable to have a theoretical and general interpreta- 
tion of how the metric on shape varies when the support is deformed. Namely, if X 
is a shape and a RKHS W on varifolds is set, we want to express the variation with 
respect to X of terms like {lJiXil^')w' with ji' a certain element of W or equiva- 
lently of juxi^) for uj G W. The proper formulation is to express the derivative of 
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/ix(<^) for small variations of X obtained by flowing a vector field from X. With 
some hypotheses on X and W, it is possible to derive a formula that generalizes 
in our context the notion first variation of a varifold studied in [2 . The result is 
summed up below : 

Theorem 4.2. Let X be an orientable compact submanifold and fix its associated 
varifold. Let v any C^ vector field with compact support defined on E and consider 
the associated one- dimensional subgroup t \-^ (j)t generated by the flow of v. Then, 
if Xf = (pt{^) 'i^s the transported manifold, we have for any C^ function (x, V) h^ 
uj{x,V) on Ex Gd{E) : 

d 



dt 



..o""*^"^ = L (^ "'"" i^) -^-1-") +/,,(-'-" + (1^1-^ 



where v^ and v~^ denote the tangential and normal part on v along X and v is the 
unit outward normal along dX. 

We refer to appendix B for the precise definition of all terms and for the complete 
proof of this formula. We will just make a few qualitative comments since there 
are a few noteworthy consequences to mention. The first important remark is that 
the variation of varifold metric is controlled only by the vector field v and not its 
derivatives. This is not straightforward since a varifold kernel K{(x^V\(x' ^V^) is 
encoding in an arbitrary non-linear way first order information through the inclusion 
of Grassmannian variables V and V . Interestingly, the result is valid even if the 
dimension or the codimension of X is different from 1. In return, we see that 
the formula involves some terms on the boundary of X, one of them expressing 
the tangential extension of X along its boundary and a second one related to the 
variation in the tangent space direction on the boundary in the normal direction 
v^ . In the interior of X, we see that the variation depends only on the orthogonal 
component v^ of v, which, in other terms, shows that the gradient of the attachment 
distance is orthogonal to the shape. 

5. Large deformation matching of unoriented shapes : an algorithm 

At this point, we have defined a theoretical background to represent and com- 
pare unoriented shapes through varifolds and kernels on varifolds. What we have 
obtained is a distance (actually a whole class of distances provided by different 
kernels) between the objects that can serve as an attachment term in practically 
any matching process. In the rest of the article, we will focus on one particular 
model of large deformations called LDDMM, which has already proved its inter- 
est for various sorts of data, from images [4 , curves and surfaces [16] [12] to fiber 
bundles [11 and more recently functional shapes [7]. The purpose of the following 
is to derive the equations needed for the numerical implementation of LDDMM on 
varifolds. In all this section, we will restrict ourselves to the usual cases of curves 
and surfaces living in the 3-dimensional euclidean space R^. 

5.1. Description of the algorithm. Let's briefly remind the principal features 
of LDDMM modelling. Given two objects 0\ (the source) and O2 (the target) of 
the ambient space E (two curves, two surfaces...), the registration problem consists 
in finding an 'optimal deformation' that transforms the source object onto the 
target one, or at least approximately. In the LDDMM framework, we consider a 
group of deformations that are diffeomorphisms given as the flow of time-varying 
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vector fields of the space. If V is a certain Hilbert space of vector fields on E and 
Ly([0, 1]) denotes the space of time- varying vector fields Vt with \/t G [0, 1], Vt G V 
and Jq Ivflydt < oo, then the group of deformations associated with V denotes Gy 
is the set of all (j)i satisfying ^g = ^d and the differential equation 

The underlying geometrical setting is to consider Gy as equipped with a right 
invariant metric induced by V playing the role of the tangent space at identity 
[24] . Matching the two objects can be formulated as finding the solution of the 
variational problem : 

(19) inf J{v),J{v) = \[ \vt\ldt^ld{^\.0^,02f 

This functional J can be interpreted as the sum of two terms, the first one that 
constrains the square length of the deformation path (which is the meaning of the 
word optimal deformation used above) while the second term drives the matching by 
measuring the distance of the deformed source (j)i.Oi to the target. 7 is a trade-off 
parameter between the two terms. On the optimality of the deformations due to the 
minimization of the first term, much has been said and proved in the past, notably 
with the interpretation as geodesies in shape spaces (cf ^). This resulted in two 
main algorithms for the resolution of the variational problem, namely a gradient 
descent scheme in the space i^y ([0, 1]) and later a geodesic shooting algorithm. Both 
these schemes are essentially related to the dynamics of the deformations and not 
to the nature of the objects and can be therefore adapted to the varifold case in a 
similar way. To fix ideas, in the following, we shall focus on the gradient descent 
algorithm for non-oriented shapes but the equations given for the data attachment 
distance and its gradient could be plugged into a geodesic shooting procedure almost 
straightforwardly. 

If the source object Oi is given as a finite set of points in the ambient space 
{q^}i=i^,,^Ni and if the vector field's space V is taken as a RKHS with kernel Ky^ 



then it has been shown in p^, that the optimal vector field in the problem (19) can 
be searched under the particular form : 

(20) vt{x) = ^i^y(^j,x)aj, 

i=l 



with ql = q^ ^ / Vs{ql)ds . 

^0 



This means that the optimal vector field at all times t is completely parametrized 
by the set of vectors a\ ^ E that we call the momenta of the deformation. The opti- 
mization of J is therefore equivalent to the minimization of the following functional 
with respect to a = {a\) : 



(21) 



^(a) = \ /o T.ii{Kv{qi, qi)4,<^\)dt + g{qi) 



This leads to finite dimensional optimal control problem where the vector at is the 
control variable. As detailed in [25], the variation of J with respect to a variation 
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5a of a gives: 






5J{a)= / '^{'^Kv{ql,4){o^l-if,),5a\)dt 

3 

where the covariable 77 is solving the backward integration scheme : 

f]t = -dv*{q)r]t 
(22) 



where dv* is the adjoint of dv and v is given by (20). Since {Ky{ql^ Qi)ij) i^ positive 
definite, t ^^ 7^^ — at as a descent direction for J. 



At this point, we must insist on the fact that the differential equations of (21) 



and (22) are focused on the dynamics of optimal diffeomorphisms, independently 
from the nature of the objects, whether they be sets of landmarks, oriented or 
unoriented curves and surfaces, measures... The only requirement is to specify the 
attachment quantity g{qi) for the computation of the total energy and its gradient 
Vg{qi) for the initialization of 77 in the backward equation. Adapting LDDMM to 
unoriented sets represented as varifolds therefore consists essentially in computing 
these terms. We refer to appendix B where we have detailed such computations for 
the case of embedded curves and surfaces in R^ compared through a kernel metric 
on the space of varifolds. The rest of the algorithm, similarly to LDDMM for 
landmarks or currents, is a gradient descent on the at computed through equation 



(22) 



5.2. Some simulations and results. We now come to a few results of varifold 
LDDMM algorithm, in which we want to emphasize the benefit in situations that 
traditionally involve orientation issues when objects are matched with the currents' 
framework. In all the following experiments, we chose the Gaussian kernel 

ke{x,y) = e ^e on the space E. As discussed previously, there are also many 
possibilities for the kernel kt on the Grassmann manifold : we have focused es- 
sentially on the Cauchy-Binet kernel and the Gaussian kernel. Even though both 
kernels were proved to induce a distance on the set of reunion of submanifolds, 
we have also explained why the Gaussian kernel has better separation properties, 
especially when the scale a^ of the kernel on E is chosen to be large. Still, the 
Cauchy-Binet kernel has the advantage of not introducing an additional scale pa- 
rameter and also leads to fast numerical computations for shapes with high number 
of points. We have therefore used the Gaussian kernel with experiments on curves 
but rather the Cauchy-Binet kernel for 3D surfaces which are usually more sampled. 
The example of figure [5] shows a situation with pikes' structures that are typically 
not well matched using LDDMM with currents. This example can be also used to 
show the robustness of varifolds when we increase the scale parameter of the spatial 
kernel /Cg, as illustrated in figure |6] It is particularly remarkable that even at scales 
much larger than the size of the object, the algorithm is still able to capture the 
overall description of the target shape, simply based on the distribution of tangent 



spaces (cf discussion at the beginning of section 4.3). 

In figure [7| we go back to the example of figure [1] where shapes are constituted of 
many disconnected pieces of curves which are not consistently oriented, and show 
the results of both methods (currents and varifolds). Other interesting situations are 
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Figure 5. Example of registration between 2D curves. The source 
shape is the blue circle, the target is the red star with four narrow 
branches. On the left, the matching is performed with the approach 
of currents. On the right, with the approach of varifolds exposed 



in section 5.1 with the same parameters. We see that the branches 
are well recovered with varifolds whereas the current's metric is 
nearly insensitive to them. 
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Figure 6. The same example for varifolds but performed at dif- 
ferent scales Ge- 
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Figure 7. Registration result on the example of figure [T] with the 
framework of currents (left) and of varifolds (right). As expected, 
varifolds can deal efficiently with any orientation of the different 
components. 

given by tree-like structures for which orientating the different branches consistently 
for registration can become a nearly intractable problem when faced with many 
intersecting branches. We give a simple insight of this phenomenon in figure [8J 
We also present some results of surface matching done with our varifold algorithm. 
The first one is a matching between two thin envelopes. Due to the proximity of 
the lower and upper membranes, current-norm based registrations have a tendency 
to squeeze both membranes together in order to eliminate unmatched parts of 
the shapes, which we can see on figure [9| Finally, in figure [Toj another surface 
registration is applied between a standard sphere and the Stanford bunny model. 
We observe again some mismatch artifacts in the currents's result : the ears are not 
fully recovered as well as some small details in the head and we see the apparition of 
undesirable membranes at the basis of the ears. In contrast, the varifold algorithm 
matches almost perfectly all parts of the shapes. 



6. Conclusion 

Having explored in the ffist place the orientation issues that come along with 
the use of currents in shape modelling, we have adapted in this paper the alterna- 
tive concept of varifolds, by representing shapes as spatially spread distributions 
of unoriented tangent spaces themselves modeled as elements of the Grassmann 
manifold. As for currents, we embed a wide variety of objects like submanifolds 
and rectifiable subsets in one common functional space. The ffist contribution of 
our work was to define metrics that share theoretical requirements and enable con- 
venient numerical computations. We have shown how this can be done by defining 
Hilbert structures on varifolds through reproducing kernels and proposed a general 
construction process for such kernels. From the theoretical side, we proved that 
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Figure 8. Example of registration on trees with currents (left) 
and varifold (right). Observe again that, with currents, consistent 
orientation between the source and the target is necessary to avoid 
the unnatural deformation shown on the figure. 



these metrics, as opposed to the previous representation by currents, do not arti- 
ficially 'eliminate' shape volume due to orientation. We also computed a variation 
formula of the metric with respect to shapes. In terms of numerical applications, 
we have presented an adaptation of large deformation shape registration based on 
this framework. The results we presented are mainly focused on synthetic examples 
for the time being, but we believe there are already important conclusions and per- 
spectives that are brought out. The first one is the fact that our approach does not 
require any orientation of the shapes, which, in our opinion, can be of great interest 
in datasets which are not easily or even not orientable in a consistent way. Another 
point we would like to suggest is the idea that the tangent space distribution of a 
shape, even if the positions are not accurate (e.g for large scale kernels on E)^ is 
already giving much information on the shape itself. Because of orientation and the 
canceling effect, currents totally loose such information. In contrast, varifold-based 
algorithms are able to recover shapes, at least approximately, even at very large 
spatial scale. 
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Source surface 



Target surface 




Matching with currents 



Matching with varifolds 



Figure 9. Registration of two heart envelopes from mouse em- 
bryos (see Le-Garrec et al, [21]). We show the source and target 
envelopes on top, and the result after matching with currents and 
varifolds. Observe that the algorithm based on currents is able to 
reduce the data attachment distance between source and target by 
crushing both sides together instead of displacing them, unlike the 
varifold framework. 



Appendix A. A variation formula for varifold metrics 



This appendix is dedicated to the proof of theorem |4.2[ We assume that X 
is a smooth compact orientable submanifold of E so that we can consider the 
canonical volume form gx on X. Given a C^ vector field ^' on £^ with compact 
support, we consider the 1-parameter group of diffeomorphisms (j)t with (pQ = Id 
and dt\t=o(l)t = v. The whole point of the theorem is to compute the variation of 
Mx(^) when X is deformed by vector field v, that is to say : 



(23) 



d 
di 



t=o 



d 



t=o ^x 



uj{x,T,Xt)dn\x) 



where Xt = (/)t(A'). Note that this generalizes the first variation of a varifold 
computed in [2^, for which cj = 1. Here, uj is any C^ function on £^ x Gd{E). As 
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Figure 10. Surface matching results between a sphere and the 
Stanford bunny. At top row, the source sphere and the target 
bunny shape. In the middle row, the result after matching with 
currents. At bottom row, with the varifold approach. 
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already explained in section l3J we have /i0^(x)(^) = Mx(0t^) and : 
(24) 



d 
dt 



t=o Jx <^^ 



t=o Jx 



where £yUJ stands for the Lie derivative of a;(x, V) along v. Hereafter, to simplify 
the notation, for a function (x, V) ^^ /(x, V), we will simply write J^ / instead of 
J^f{x,T,X)dn^{x). Now, {rtOj){x,T,X) = \d,c^t.T,X\uj{U^),d,MT,X)). Let 
us recall that Jt = \dx(t)t-TxX\ is the volume change in the direction of the tangent 
space TxX. Taking the derivative inside the integral at t = leads to three terms : 
differentiate the function uo with respect to position, with respect to the tangent 
space direction and differentiate the volume change. The first term is the simplest 

one and an immediate computation shows that it equals ( ^— |^ ) • The two others 

are more involved. 

Derivative of the volume change : For any vector field u defined on X, we 
shall denote by v7 and u^ the tangential and normal components of u with respect 
to the tangent space of X at each point. We also introduce the connection V.- on 
the ambient space and an orthonormal frame of tangent vector fields (6^)^=1^. .^^^ on 
X. Now Jt = Y^det(((ia;(/)t(ei), dx(l)t{ej))ij so a simple calculation shows that : 



dt 



t=0 i^i 



Writing v = v^ -\- v-^ provides a first term X]i=i(^i5 ^e^'^^) which is the tangen- 
tial divergence of the vector field v^ denoted usually divx('^^). The second term 
becomes X]i=i(^*5 ^e^'^^)- For all i = l,..,(i, we have (ci^v^) = so that after 
differentiation we find that (e^, Ve^'^^) = — (Ve^e^, i'^). Therefore : 

d d 




\i^l / / 

In this last expression, we recognize the mean curvature vector to the subman- 
ifold X, which is the trace of the Weingarten map and is denoted Hx- As a result, 
we find that : 

/a;— Jt= ujdivxiv^)- / uj{Hx,v^) 
Jx "^ ^=0 Jx Jx 

Now, we will show that the first term can be rewritten as a boundary integral. 

Indeed, if we denote by Q the function defined on X by uj{x) = uo^x^T^X)^ we 

have divx(^'^^) = a;divx('^^) + (Va;|'U^). Applying the divergence theorem on 

the orientable manifold X gives : 

/ a;divx(^^) = - / (VcD|^^) + / uj{u,v~') 

Jx Jx JdX 

where v is the unit outward normal to the boundary. 
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Derivative with respect to tangent spaces : We now come to the last term 



in equation ( 24 ) , which is the variation of uj with respect to the tangent space in 



Gd{E). As explained in the beginning of section 3.1, we can identify the tangent 
space to Gd{E) at V with the space of linear applications C{V^V^). In addition, 

if we set Vt = dx(j)t{TxX)^ we have : 

d 



dt 



Vt = Pt.x^ o \/v\t^x e C{TxX, (T,X)^) 



t=o 



which we will note more concisely 



dt 



Vt = V^v. The space /:(T^X, (T^X)^) 



^=0 



being trivially isomorphic to T^X-^^iT^X)*^ we can consider V^v being an element 

duj 
of that space and we introduce 7-— which we therefore identify to an element of 

oV 

{TxX^Y (g) TxX : — - = YJj=d+i Vj ^ ^j foi" {Vd+i^ "^ Vn) an orthonormal frame of 

TxX^ and (a^) vectors of T^X (as usual r]* denotes the linear form (77, .)). Then 
the variation we wish to compute is : 



dV 






duj , 



If we introduce ( ^|^ ] = Yl'j=d+i ^jM^j = YJj=d+i{^3->^)^J which is a tangent 
vector field on X, we have : 

divxf^l^j =Y^ Y^ {{ei,Ve^oij){r]j,v) ^ {ei,{Ve^r]j,v)aj) ^ {ei,{r]j,Ve^v)aj)) 

^ ^ i=l j=d-\-l 

The last term in the sum is also X]7=d+i^^i' ^"j^)' which is nothing else than 



-—-\\/v ) . As for the two other terms in the sum, it's easy to see that it equals : 



^ d n \ . 

^(e„Ve, Yl ^I^^^-)l^ =(divx( 



dV 



\v) 



So get eventually that : 
duo , 



(25) 



^|V.l=divxl^l. 



duj , 



(divx ( ^ ) 1^) 



Integrating equation p5|) over the submanifold X and using as previously the di 



vergence theorem, we find that : 



(26) 



\v) 



Synthesis : Grouping all the different terms obtained so far, we get the follow- 



mg : 



I/"" = X(^-^'-(^)i'')-X<^*^>-"<«-i""> 



' duj 
IX \9x 



- UJV 
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We remind that Cj{x) = uj{x^TxX) so {\/Cj\v^) = I -^|'^^ ) + I -^tt|V'U^ ) and 



applying the result of equation (25) to v^ : 

duj 
dV 



V.-)=div,(|i|.-)-(div.(^)|,-) 



Using the divergence theorem on div^ ( ^77 1 ^^ ) ^^^ the equality v = v~^ -\- v-^ we 
find eventually that : 



which proves the result of theorem 4.2 



Appendix B. Discrete computations for varifold LDDMM 



As explained in section 5.1, the LDDMM algorithm relies on the computation 
of a distance between the deformed source object ^(C^i) and the target object (^2, 
and on the gradient of this distance with respect to final points qi of the object Oi. 
In our context of unoriented shapes modeled as varifolds, we use distances provided 
by kernels and their associated RKHS presented in section |4j Those are the tensor 
product of a kernel ke on the space R^ and a kernel kt on the Grassmann manifold. 
If W is the RKHS corresponding to k = ke ^ kt^ the attachment distance we have 
is then : 

(27) d(0^Ol,O2) = l|(<^D*m -MoJl'v^. = W^^rao,) - f^o.fw^ 

Since ^^(Oi) is the transported source object by the vector field at time 1 (repre- 
sented by the set of points qi with the previous notations), formally the problem 
reduces to express explicitly quantities like ||/ix — /iy||^* where X and Y are two 
shapes of the same dimension and compute the gradient of such terms with respect 
to the points of X. 

We assume now that X and Y are two unoriented (i-dimensional shapes (possibly 
disconnected) given respectively as sets of vertex {xk)k=i,..,N^ {yi)i=i,..,M together 
with sets of unoriented simplexes (//, .., //)i=i,..,n and (^j, ..,^j)j=i,..,m- The com- 
putation of the varifold representations /ix and /iy of each shape was explained 
in section [3^1] and provided by equations ^ for curves, (IgI) for surfaces. We write 
fix = '^i=ik-^{pi,Ui) and /iy = Xl^i Aj.(5(g^. y^.). Then the attachment distance 
becomes : 

^ = IImxII^* -2(/ix,/iy)H^* + ||/iy||^* 
= ^ liljke{pi,Pj).kt{Ui,Uj) 

i,j=l..n 

-2 ^ ^ li\jke{pi,qj).kt(Ui,Vj) 

i=l..n jf = 1..77T, 

(28) + Y. ^i^jke{qi,qj)MVi,Vj) 

i,jf = l..n 

Thus, computing the distance between X and Y consists basically in computing 
the varifold representation of the shapes from their sets of points and meshes and 
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then make repeated kernel evaluations. If ke and kt are specified, the previous 

expression becomes totally explicit. 

The more technical part is the computation of the gradient of the previous distance 



with respect to the (x/c)'s, the points of the first shape. In equation (28), the 



attachment is a function of the p^, li and Ui but all these terms are themselves 
functions of the (x/e)'s, respectively as the centers, lengths and tangent directions 
to every cell of X. The resulting function that we have denoted g{{xk)) can be 
therefore differentiated as a composition : 



(29) 



9xk9 = yZi^P^^^^XkP^ 



i=l 



• di^A o dxji + du,A o dxj^Ui) 



The differentials of A with respect to the p^, k and Ui are easily expressed from 



equation (28), involving the differentials of the kernels ke and kt. 

' 9p^^ = Yl]=ihlj[dike{pi,Pj) ^ d2ke{pj,Pi)].kt{Ui,Uj) 
-2 YlJLi kXjdikeiPi, qj).kt{Ui, Vj) 



(30) 



9u,A = YJj^i liljke{pi,Pj).[dikt{Ui, Uj) + d2ke{Uj, Ui)] 
-2 Y.'jLi li><jke{Pi, qj).dikt{Ui, Vj) 



duA = 2YTj^iljke{PuPj)'hiUi,Uj) 
-2 E^i >^jke{Pi, qj).kt{Ui, Vj) 
where di and 82 denotes the derivative with respect to the first and second argu- 



ment. We remind (cf section 4.1) that, in our approach, kernel kt is the restriction 



of a kernel defined on the vector space C{E) so that writing its derivatives is nat- 
ural. 

Now, the last step consists in computing derivatives of p^, li and Ui with respect 
to the points of X. These depend on the dimension and codimension of the shape. 
In the case of curves and surfaces embedded in R^, it can be simply done by differ- 
entiating equations (|5| for curves or equations (|6| for surfaces. We give the details 
for these two cases below. 

' \xp — xji|. The tangent space 



In the case of curves, we have pi 



H — \^ 



direction can be represented as a single normalized vector Ui = 
gives for all i G {1, .., n} and k G {1, .., N} 



•' i •' i 



which 



(31) 



dx,,Ui 



1 



Xf2 



J 1 ) 



Xfi 



fh^ 



3 /• 



\Xf2 


-Xfi\ 


1 


/ 


l^/f ~ 


^//l V 


^ff^S 


~^fl^s 



\Xp — Xfi\ 



Hk=ff} 




{k=m 



\Xf2—Xfi\ ' -"- ' \Xf2—Xfl\ 

In the previous equations, x^^s denotes the number s coordinate of x^ in the em- 
bedding space's canonical basis (e^). These equations can be simply interpreted 



32 NICOLAS CHARON, ALAIN TROUVE 

and implemented as the way to distribute the gradient computed with respect to 
the segments over the points of the shape X. 

In the case of triangulated surfaces in Mr ^ we have pi = — ^ ^ ^, 

h = 1(^/2 — Xfi) A {xf3 — Xfi)\ and again the tangent space is simply encoded by 

the normalized no 
and k G {!,.., A/"} 



{x.i-Xj,2)A(Xr2-Xj,3) 

the normalized normal vector m = r? — ^ Hr? — ^ ^ so that for all i G { 1 , . . , n | 

* \{X ^l-X ^2)/\{x ^2-X ^3)\ I •> •> } 



(32) 



dx^,sPi = o (^{/c=//} + hk=m + hk=ff}) 



dx^.. 



^Ui= - Us /\ {Xf2 - Xf3) - Us A {Xf2 - Xfs) , Ui\ uA ^{/c=/,i} 



dxkji = (^s A {Xf2 - Xfs) , ii^^ ^{k=fl} + ••• 

We wrote "+..." to mean that similar terms are obtained for k = // and k = ff. 
Combining the last relations with equations ( 29 ) and ( 30 ) provides a full description 
of the gradient computation in practice. 
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